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Abstract 



We study equilibrium properties of a system of particles in two dimensions, 
interacting with pair and three body potentials, which undergoes a structural 
transition from a square to a rhombic lattice and thus constitutes a simple 
model for a generic tetragonal to orthorhombic transition. We aim at an in- 
termediate level of description lying in-between that of coarse grained elastic 
strain hamiltonians and microscopic ab -initio approaches. Macroscopic ther- 
modynamic properties and phase diagram at zero and finite temperatures as 
a function of the density and the relative strengths of the pair and three body 
energies are obtained using lattice sums, an approximate 'cell-model' theory 
and molecular dynamics simulations in the NVT ensemble. We propose that 
this model solid can be used as a test bed for studies of statics and dynamics 
of structural transitions. 



I. INTRODUCTION 

In spite of its fundamental and technological interest, there is as yet no general theoreti- 
cal framework to predict the final microstructure of a solid following changes in temperature 
or stress across a structural transition IQ. This is in part because most experimental stud- 
ies have focussed on technologically important solids which are far-from-idealized ; thus it 
has been difficult to isolate generic principles amidst the volume of empirical data P,^- 
Furthermore, the spatial and temporal resolution of in-situ experimental probes is limited ; 
this makes it difficult to follow microstructural changes at short length and time scales. We 
believe progress can only be made if we (a) identify simple 'model' systems which would 
serve as an arena for detailed studies on the dynamics of structural transitions in solids 
and (b) develop 'probes' to study the dynamics and morphology changes at short scales ; in 
other words to follow the motion of individual atoms as the transformation proceeds. High 
speed computational modelling allows one to make useful contributions to both (a) and (b). 
This paper concerns point (a) - we provide a model system for a simple structural transition 
which could serve as a test bed for future studies of statics and dynamics of structural tran- 
sitions. In another paper 0] we have discussed point (b) - we have studied the dynamics 
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of solid state transformations in this model system by tracking individual particles as the 
transformation proceeds. 

Our choice of the simple structural transformation, a square to a rhombic lattice in two 
dimensions (2-d), is motivated by two important considerations. The first is its relevance 
to 'real' solid state transformations. The square-to-rhombus transition may be regarded as 
a rather accurate representation of the three dimensional tetragonal to orthorhombic 
(TO) transition |^ in an oriented single crystal where the strain along the third direction (c 
- axis) is negligible. This transition is observed in many technologically important systems 
1^, such as the high-T^ compound YBa2Cu307. 

The second motivation is more conceptual and warrants some explanation. Following 
a quench across a structural transition, the atoms constituting the solid have to rearrange 
themselves, since the parent phase is thermodynamically unstable. Instead of moving to the 
new equilibrium configuration, the motion of the atoms is arrested in a final microstructure 
which is very far from equilibrium. The microstructures often display features at length 
scales ranging from lOOOA — lOO/im, many orders larger than the lattice spacing. A more 
appropriate description of the dynamics at these scales is in terms of continuum degrees of 
freedom. It is however not clear apriori, what are the relevant continuum degrees of freedom, 
especially in situations where the solid undergoes large deformations from the parent. In 
most theoretical studies of this problem, the only continuum degrees of freedom relevant at 
these scales have been taken to be the components of the strain tensor HJToD . However, short 
length scale phenomena like atomic rearrangements flllJl^ are not captured by these strain- 
only theories and may in some cases affect the kinetics of the transformation. The only 
unbiased way to determine the complete set of relevant degrees of freedom in the continuum 
is to start from a more microscopic description and arrive at a continuum description using 
a coarse-graining method [|T3|. Such a coarse-graining program is more easily set up in the 



crystallographically simpler square-to-rhombus transition. 

How microscopically detailed should our microscopic model solid be ? An ab-initio or 
semi empirical description [0] which includes electronic degrees of freedom tailor-made for 
a real system such as YBa2Cu307 suffers from 3 drawbacks — (i) it compromises the need 
for generality, (ii) it is computationally expensive and (iii) it is difficult to extrapolate to 
the continuum. We therefore model an effective Hamiltonian accurate over distances smaller 
than the bulk elastic correlation length but larger than the typical atomic spacing. This 
effective Hamiltonian, coarse-grained over the faster electronic degrees of freedom, will in 
general have pair and many-body interactions. 

Our model system therefore consists of a set of N 'particles' confined in a 2-d box of 
volume V at a fixed temperature T. The particles interact with each other via specific two 
and three -body interactions. While the 2-body interaction stabilizes a rhombic lattice 
(triangular for isotropic two-body potentials), the 3-body interactions have been constructed 
to favour a square lattice. Unlike in a real solid, the square to triangular transition in our 
model solid may be driven by independently (i) increasing the density, (ii) decreasing the 
temperature or (iii) decreasing the relative magnitude of the three-body potential. We 
confine ourselves here to the equilibrium aspects of this transition and postpone a detailed 
study of the nucleation and coarsening dynamics to a subsequent publication 

In the next section we introduce our model, describing the pair and three -body poten- 
tials. In Section III we present the static lattice (zero temperature) results for the energy, 
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stress and the elastic constants and exhibit the zero temperature phase diagram. In Section 
IV we discuss the effect of finite temperatures. We present results both from a "cell model" 
approximation p!6| , p!7[] and from molecular dynamics simulations [|I8| in the NVT ensemble. 
In Section V we discuss a simple generalization of the model to include a molecular solid 
with a complex basis motif. Section VI presents a summary and conclusion of this work. 



II. THE MODEL POTENTIAL 



In this section we motivate the form of the effective hamiltonian, coarse-grained over the 
faster electronic degrees of freedom. Instead of starting out ab-initio, we shall assume that 
the solid can be described by a general functional of the densities of particles. A simple form 
of the density functional, proposed by Ramakrishnan and Yussouff (RY) |^ , views the solid 
as an extremely inhomogeneous liquid with a non-uniform, periodic, coarse-grained density 
p(r). The advantage of this approach is that it allows us to make accurate statements at 
any finite temperature ksT = I3~^. Following RY, we may write the Helmholtz free energy 
Fs per unit volume V of the solid as. 
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The leading term gives the ideal gas contribution (the integral is over the Wigner-Seitz unit 
cell of volume v), while the subsequent terms arise from the interactions between density 
waves pg with wave- vector G belonging to the set {G} of reciprocal lattice vectors (RLV) 
of the crystalline solid. Note that we have defined the Fourier transform as. 



Pg = -[ d^re 

V Jv 



iGY 



p(r). 



(2) 



The interaction terms involve the Fourier transform of the n^'^-order direct correlation func- 
tions C*^"^(Gi, G2, ■ ■ ■ , G„_i). These functions may either be evaluated from a liquid state 
theory ||22| or deduced from scattering experiments (for instance, C^'^\q) = 1 — pq=o/S{q), 
where S{q) is the structure factor of the liquid at wavenumber q) ; however evaluation (or 
measurement) of direct correlation functions beyond the second order is extremely diffi- 
cult |]23[|. Approximate ways of incorporating the effects of these higher order correlations 
have been used with varying degrees of success ||2^, though often a simple weak coupling 
(mean-field) approximation, C'-"^ = V^'^'^ /ksT, where y^""-* is the n— body potential, works 



remarkably well 25 



While the ideal gas term always prefers a uniform liquid, the sign of the interaction 
term decides whether the interactions stabilize (destabilize) density waves with wave-vectors 
G 7^ 0. As an example, let us consider first the effect of only the second order terms in 
Eq.|l|. The coefficients, C(2)(|G|), measure the stability of a density wave with wavenumber 
|G|. In Fig. |VIII| we have plotted C^'^\q) against q for a slightly supercooled hard disk liquid 
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26| in two- dimensions. This function is oscillatory and has a primary peak at roughly the 
wave number corresponding to the magnitude of the smallest RLV of the thermodynamically 
stable solid. Choosing the lattice parameter such that the smallest RLV coincides with the 
first peak of C*^^^(g), we have plotted the positions of the RLVs for the triangular (top) 
and square (bottom) lattices in Fig. |V111] . Closed packed lattices have RLVs which are, on 
an average, more widely separated than those of an open lattice. As a consequence, open 
lattices often have RLVs lying in the region of the first minimum of C'-^^(g) which is negative, 
thereby contributing to a destabilisation of the lattice. From Fig. |VIIl] we observe that this 
is indeed the case for the square lattice in two-dimensions Density waves corresponding 
to the second RLV of the square lattice with Miller index {11} are not favoured, making the 
square lattice unstable in two-dimensions. 

Note that the discussion above is the finite T generalisation of the zero temperature 
result that for isotropic, purely repulsive pair potentials in two dimensions, one can only 
stabilize the close-packed triangular lattice [|T^. For instance, a static calculation of the 
T = elastic moduli of the square lattice reveals a shear instability which spontaneously 
distorts the lattice till it regains a close-packed structure. 

In order to stabilize the square lattice one needs to go beyond the second order contribu- 
tion and consider the effect of three-body correlations. As is clear from Eq. |I], a positive (and 
large enough !) contribution from C^^\Gi, G2, — Gi — G2)p(Gi)p(G2)p(— Gi — G2), where 
any one (or two) of the wave- vectors equals {11}, can compensate for the destabilizing effect 



of the second order correlator [^]. There are many choices for the wave- vectors G, but the 
simplest combination is Gi = G2 = {10} (so that {10} -|- {01} = {11}). A straightforward 
way to ensure that this combination of density waves is stabilized is to stabilize the real space 
triangle involving the direct lattice vectors (10), (01) and (11) (and those related to them by 
symmetry). Within a simple minded mean-field approximation this may be accomplished, 
as shown below, by choosing an appropriate three-body potential which favours 0°,45° and 
90° bonds. Higher- order interactions involving four or more particles, though present in 
principle, are not necessary for our purpose. 

Our model system, constructed from this level of coarse-graining, therefore consists of 
'point-particles' interacting with effective pair and three -body potentials. The interaction 
energy E of the system is given by, 

^ = ^E^2(r.,) + i J2 ^3{rij,rjk,rki) . (3) 
For the pair potential we take, 

= V,{-^)'\ (4) 

which is purely repulsive and therefore the system has to be confined with a uniform hy- 
drostatic pressure (see below). A purely repulsive system simplifies our analysis since there 
is one fewer length scale and one fewer non -solid phase. Without the 3-body potential, 
our system is characterized by only one parameter instead of the two (temperature and 
density). On including the 3-body potential, we lose this simplification, but the variation 
of thermodynamic properties with density is still weak. Without loss of generality we can 
take V2 and a to be our units for energy and distance respectively. 
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The 3-body potential is 




+ fjk sin^(4efc)/fci 



(5) 



where the function, 



fij = f{nj) = {rij - r^f 
= 



Tij < To 

otherwise 



(6) 



and we have used the notation = |rjj|. The angles are as defined in Fig. 2. The function 
fij provides a cutoff for the 3-body potential ; as long as fij is short ranged, the actual form 
of this function does not affect the qualitative results. 

It may appear that a three-body potential requires a large investment in terms of com- 
puter times. This apprehension is fortunately unfounded. The form of this potential ensures 
that three-body energies can be calculated |T5[ extremely efficiently, requiring a computa- 
tional effort not exceeding that for the pair part. This is discussed in the Appendix. 



At zero temperature, the equilibrium configuration is a solid which minimises the po- 
tential E. Since we work in the constant NVT (and shape) ensemble, the density p = N/V 
is a constant. Assuming that the only minima of E correspond to the triangular or square 
phases, we have numerically deduced the T = phase diagram in the p — V3 plane (Fig. 
Later in this section we show, within a restricted variational calculation, that these are the 
only minimisers of E. As we see from Fig. ^, the triangular lattice is the lowest energy 
phase at high densities and small values of V3, i.e., wherever the pair interaction dominates 
over the three-body part. Across the boundary there is a strong first-order transition with 
a discontinuous change in the slope {dE / dp)y,y 

To deduce the nature of the order parameter distinguishing the square from the triangular 
phase, we look at how a rhombic lattice may be obtained from a square. Such an analysis 
makes contact with continuum elasticity in a natural way. 

At zero temperature, a continuous family of perfect rhombic lattices (labelled by posi- 
tion vectors R"^) can be obtained from the perfect square lattice (labelled by R*^) by the 
transformation, R"^ = (I + T)R° where the transformation matrix T is. 



The parameters {a = 1,2,3) are related to the components of a strain tensor by the 
following construction. We choose to measure all distortions and energies with respect 
to the undistorted square phase — our reference state. The microscopic displacements 
Ur,o = R^ — R" are therefore defined at every HP on the reference lattice, i.e., we use 



III. ZERO TEMPERATURE RESULTS 




(7) 



Lagrangian [0,^ coordinates. The full nonlinear Lagrangian strain tensor e^j is. 
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^ fduj ^ duj ^ duk duk \ 
2 \ dvj dr-i dr-i dvj J 

where the indices i,j go over x and y. The parameters ea in Eq. |^ represent the combinations 
^xx + ^yy, ^xx ~ ^yy ^-nd e^y respectively, which reduce to the usual volumetric, deviatoric and 
shear strains once nonlinearities are neglected. Note that one may start with a prescribed 
square and end with a final triangular lattice in more than one way — the transformation 
parameters are not unique. For instance, for a given orientation of the parent square lattice 
shown in Fig. Q^i) one obtains a rhombic lattice using 62 = and tan 6* = £3/(1 + ei). 
Equivalently (Fig. 1(11)) the square lattice (first rotated by 45°) can be transformed to a 
centered rectangular lattice with €3 = and 6/a = (1 + 62 + ei/2)/(l — €2 + ei/2). The two 
transformations are completely equivalent. 

One of the ofshoots of this nonuniqueness, is that any rhombus obtained as a uniform 
deformation of a perfect square, can be represented by two independent parameters ei and 
€3. In addition, since the density p is a constant in our NVT ensemble, ei is related to 63; 
ei = ^1 + 63 — 1 ^ £3/2. Thus all rhombic lattices considered by us can be labelled by a 
single parameter €3 = e (which by definition is for the perfect square lattice). This makes 
the shear strain e a good order parameter which distinguishes the square from the triangular 
lattices. 

We may now calculate the energy of T = configurations as a function of the order 
parameter e. A calculation of the energy E and its derivatives (elastic moduli) for a given 
lattice involves computation of lattice sums. We start with a finite square lattice containing 
10 X 10 sites and allow the transformation T to produce a continuous sequence of rhombic 
lattices labelled by e. We have checked for convergence of the lattice sums by increasing 
the size of the lattice and observing the consequent change in the numerical values. In 
Fig.^ we have plotted the energy per particle E/N as a function of the parameter e for 
various values of V3 (keeping p fixed). Note that for large values of V3 there is only one 
minimum at e = so that the square lattice is the only stable phase. For smaller values 
of V3 two additional degenerate minima appear at e = ±eo = ±0.27812 which correspond 
to the triangular lattice. The transition is first-order with the order parameter jumping 
disco ntinuously, |Ae| = eo. Note that in this model the change in the shear strain across the 
transition is fixed ; we shall return to this point in Sect. V where we propose a variant of this 
model in which the jump in the shear strain across the structural transition can be made 
arbitrarily small. As expected, the square and triangular lattices are the only minimisers of 
E within this parametrisation scheme. 

To make contact with elasticity theory we may compute stresses and elastic moduli, 
obtained by evaluating appropriate derivatives of the energy keeping T, N, V constant — 
dE/dea = ctq, (Fig. ^, and d'^E/deadep = Cap (Fig.0). Note that our system is always 
under a hydrostatic pressure P = (Xi ; the constant density constraint implies that for e 7^ 
there is an applied shear stress 0-3 = P(e)e (Fig.|^). This implies that the slope of the shear 
stress vs £3 curve is not the shear modulus C33 (defined for zero external stress) but C33 + P 
(Fig. 0). 

At this stage, we find it useful to point out that the results of Figs, ^^can be rationalized 
using a systematic power series expansion of the energy in terms of e^. Although such 
expansions are quite common in the literature our results show that fourth order 

terms in e^, especially cross couplings of the form e\e\ and e\e\, together with coupling to 
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the external hydrostatic pressure P need to be included in order to reproduce the T = 
results accurately. The coefficients of all these terms are however not independent. For 
instance, relationships like, 



_ ^ 
eo del 



(9) 



dictated by the geometry of the triangular phase have to be satisfied for all temperatures 
and densities. 

In this section we have been able to show that our model potential indeed produces 
the square and the triangular lattices as minima of the energy. The potential parameters 
maybe tuned, if necessary, to a real system by comparing the elastic properties of this model 
system to experimentally measured quantities. By varying the density or the strength of 
the three-body potential one obtains a zero temperature first order structural transition 
between a square and triangular lattice. What happens to the structural transition at non- 
zero temperatures? We study this question in the next section. 

IV. NONZERO TEMPERATURE RESULTS 

In this section we analyse the phase diagram at T 7^ as a function of V3 or p. We do 
this by two methods — an 'exact' molecular dynamics (MD) simulation in the constant 



NVT ensemble, using the 2- and 3-body potentials defined earlier, and an approximate 
'cell-model' fl^ based on the deformation parameter e. The latter leads naturally to an 



approximate continuum elasticity description at T 7^ 0. We take up the cell-model analysis 
first and compare its results with the exact MD simulation in the next subsection. 

A. Cell model approximation: free energies and phase stability 

Imagine being in a region of the zero temperature parameter space V3 — p, where the 
square solid is the stable minimum of the energy. As the temperature is gradually increased, 
the contribution of the phonon entropy to the (Helmholtz) free energy destabilizes the square 
lattice. In order to quantify this effect one needs to go beyond the static lattice and con- 
sider phonon fluctuations. Although a direct calculation of the contribution of phonons to 
the lattice energy is straightforward pO[, we choose to use the much simpler, though not 



necessarily less accurate 'cell-model' approximation. 

Before discussing the cell- model approximation, let us mark its regime of validity. First, 
the cell-model approximation neglects contributions from topological defects like dislocation- 
antidislocation pairs and thus breaks down near the melting point In two dimensions. 



there is a further complication, since fluctuations of the displacement field u due to phonons 



diverge logarithmically |]19[ with system size. This divergence is however weak and may be 
ignored for the system sizes under consideration. 

Recall that at T = 0, the configurations of the perfect rhombi were parametrised by a 
single deformation variable e. Is this true at T 7^ 0, when the lattices are not perfect due to 
phonon fluctuations ? It turns out that the constraints of rhombic symmetry and constant 
density still allow a parametrization of the T 7^ configurations by a single function e(x). 
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at least when the temperatures are low. Thus the energy E may be written as a functional 
of e, which at low temperatures may be replaced by its mean (e). 

Within the cell - model approximation, the partition function of a lattice of particles 



at temperature T is given by ||T6| , 



Z{{e),T,N) = [A-' f drexp{-5U^)/kBT)f x 

exp{-E{{e))/kBT) (10) 

where A is the thermal wavelength and S(f) is the change in potential energy of a single 
particle as it moves around within a unit cell of size in a potential well arising from its 
interaction with all its neighbours. A further harmonic approximation for 6(j) leads to the 
familiar Einstein approximation. At the other extreme, for the hard disk potential, 6(f) = 
except where overlaps occur and the cell-model approximation becomes identical to the free 
volume Jl^ theory. The Helmholtz free energy for any rhombic lattice labelled by (e) may 



now be obtained by using F = —ksTlogZ. 

Evaluation of the Helmholtz free energy (Fig. H) allows us to calculate the V3 — T phase 
diagram at any density p as also the limits of metastability of the square lattice. 



B. Molecular dynamics simulations 

To obtain accurate results for the phase stability at non-zero temperatures, we have used 
a molecular dynamics simulation for our model system. We simulate A^ = 2499 particles 
(50 X 50 unit cells with vacancy to improve the kinetics) in the NVT ensemble using a stan- 
dard Nose-Hoover thermostat. Starting from an ideal square lattice, we have equilibriated 
systems at various values of V3 and temperature for a fixed density for about 50 — 100 x 10^ 
molecular dynamics time steps or till thermodynamic quantities like the pressure and energy 
have stabilized. The final structure is then examined and this information used to obtain 
the phase diagrams shown in Fig.|^. We display the phase diagram for two densities p = 1.05 
and 1.1. Together with the molecular dynamics results we have also plotted the results of 
the cell-model approximation. We observe that for low temperatures, the cell-model ap- 
proximation faithfully reproduces the actual phase boundary while at higher temperatures 
it begins to deviate. The cell-model approximation is also used to plot the limit of stability 
of the square phase in the triangular region. 

Both the molecular dynamics simulations and the cell model calculations predict that 
the square to triangular transition remains first-order over a wide region of parameter space 
even at non-zero temperatures. For larger density, the transition point shifts to higher values 
of V3. This is expected since a high density favours the triangular lattice. Also, the square 
phase becomes unstable for lower values of V3 as the density is increased. The jump in the 
order parameter remains fixed at |Ae| = cq all along the transition line. This aspect of 
our model is specifically addressed in the next section, where we show that inclusion of an 
anisotropic pair interaction allows one to tune the order parameter jump all the way to zero. 

We end this section with the following observations. We have seen that the exact MD 
and approximate cell model give qualitatively similar results. More sophisticated phonon 
fluctuation calculations may even produce quantitative aggreement. Figure || suggests that 
the Helmholtz free energy may be expanded in powers of e, just as was noted at T = 0. 
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Though we do not exphcitly demonstrate here, we may recover elasticity theory (including 
corrections arising from thermal fluctuations) by constructively coarse-graining as in Ref. 



V. GENERALIZATION TO MOLECULAR SOLIDS 

The model 2-d solid discussed in the preceding sections has the virtue that it is simple 
enough to begin a detailed theoretical study of both the equilibrium and dynamical features 
of the TO transition across a range of length and time scales. However if we were to compare 
the results of such a study with experiments on realistic systems, we would immediately face 
a problem. Most sohds undergoing a TO transition |^,||, for instance YBa2Cu307, have a 
complex basis, consisting of many atoms per unit cell. These systems generically have much 
smaller jumps in the shear strain at the TO transition compared to the jump computed 
in the previous section. To appreciate the quantitative discrepancy, recall the discussion 
following Eq. H, where we showed that any perfect rhombus obtained as a deformation of 
a square may be parametrized by either e or b/a. Deflning an orthorhombic distortion as 
D = (b — a)/{b + a), we find that D = (\/3 — 1)/ (VS + 1) ~ 0.27 for our model square-to- 
triangle structural transition - significantly larger than D = 0.0085 for the TO transition in 
YBasCuaOy. 

Apart from this, there might be a more serious qualitative mismatch between our model 
solid and real systems undergoing a TO transition. Changes in temperature or pressure 
(hydrostatic or chemical) may lead to a local structural rearrangement (optical modes) 
which would couple to the strain tensor. The jump in the shear strain across the structural 
transition may therefore, unlike in our model solid, vary along the phase boundary, even 
going to zero (phonon softening) at a critical point . 

We shall see that we may address both these issues within an anisotropic variant of our 
model solid. Our attempt will be to incorporate the complex basis, with many atoms per 
unit cell, into an effective hamiltonian between 'point particles'. In the spirit of an effective 
hamiltonian, we will coarse-grain the density over a length and time scale corresponding 
to the 'size' ^ and relaxation time r of the basis. Thus we may define a coarse-grained 
density as p(r) = Y^fj, Pfj,{r), where /i = labels the atomic species making up 

the basis. This coarse-grained density profile p(r) will have peaks at the centre of mass 
of each basis, falling off to zero over a length scale ^ and having a cross section which is 
spatially anisotropic. If we assume that this anisotropic cross section has a fixed shape at 
a given temperature and pressure (true when the associated optical branch is much higher 
than the acoustic branches), then we may write the effective hamiltonian as arising from 
a collection of 'point particles' interacting via an anisotropic potential. The form of the 



effective hamiltonian may also be motivated in terms of a density wave picture |^0[ in much 
the same way as in Section II. An anisotropic density interacts via an anisotropic direct 
correlation function. Within a mean field approach this reduces to a pair potential which 
depends not only on the distance between the two basis motifs but also on their orientation 
relative to the crystal axes — orientation fiuctuations within the motifs being neglected. 

For specificity if we assume 2 atoms per basis, we may then arrive at the following 
modification of the two-body potential using the arguments outlined above. 
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*2(|r.,|) = V^2(^)'' X (l + asmH{e,,-^)) , (11) 



where the anisotropy parameter a has a fixed value at constant T and P. On the other hand, 
the 3-body potential may be taken to be the same as in Eq. |^. Setting 7 = 1 + a, we see that 
7 is always positive. All angles are measured with respect to the [01] axes of the undistorted 
square lattice. The angle ip represents the orientation of the basis and 9ij = sm~^{xij/\rij\), 



see Fig. |T0|. The total energy is a function of ip so that uniform rotations of the basis with 
respect to the crystal axes cost energy (optical mode) while simultaneous rotations of the 
basis together with the crystal axes is a symmetry of the hamiltonian. 



Using the modified two-body potential (Eq.[TT|), we compute the energy as a function of 
e as in Section III. For a fixed a, the total energy minimized with respect to ip and e leads 
to ip = (independent of e). As before there are 3 minima in e, one at e = (corresponding 
to a square) and the other two corresponding to rhombi with |e| being smaller than the 



value for the perfect triangular lattice ~ .28 (Fig.jTOl). The jump in the value of the shear 
strain e across the structural transition is therefore smaller than that obtained in Section III. 
Moreover we find that this jump in e goes to zero and the region over which the square phase 
is metastable shrinks and disappears as V3 ^ thus indicating a continuous transition at 
a tricritical point. One expects, therefore, that for real systems fiuctuation effects near the 
T-0 transition would be more pronounced. This fact is actually bourne out by experiments 



The zero and finite temperature phase diagrams are shown in Fig. |Tl|. The zero temper- 
ature phase diagram clearly shows the location of the tricritical point where the jump in 
the order parameter vanishes. The effect of finite temperatures is addressed easily within 
a cell-model approximation. The calculation may be carried out along the lines outlined 
in the last section. Once again we see that the Helmholtz free energy can be written as 
a power series expansion in e. The results of the calculation are expected to be accurate 
at low temperatures if the anisotropy is not too large. For larger anisotropics the effect of 
(tri)-critical fluctuations may alter the results of our simple mean-field estimates. Our result 
shows that the square lattice now becomes stable over a much larger range of V3 than in the 
isotropic case. The region of metastability of the square lattice however decreases and the 
first-order transition is weakened. 



VI. SUMMARY AND CONCLUSION 

In this paper we have described a model system which is designed to undergo a square to 
rhombus transition in two dimensions. We believe that our study will be useful in two ways. 
On the one hand, it may be used as a simple simulational model for the T-0 transition in 
real materials which often consists of a large number of individual atomic species making 
it difficult to study using ab-initio methods. For this purpose, the parameters V3 and a 
have to be "fitted" to observed properties of the particular realistic system. On the other 
hand, we could use this system to study, in general, the dynamical pathways of a simple 
first order solid state phase transition involving a structural transition. It is this context 
that we would like to emphasize. Once the equilibrium properties are determined, we look 
at the nucleation dynamics, growth modes and microstructure of the rhombic phase growing 
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in the matrix of the parent square lattice [Q. The effect of defects such as vacancies and 
dislocations are automatically incorporated in our microscopic approach. In future, we hope 
to obtain atomistically detailed information about the statics and dynamics of solid-solid 
interfaces. 
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VIII. APPENDIX 



In general, evaluation of three body energies requires sums over all possible triplets which 
for a system of reasonable size is prohibitively expensive. The particular form for the three 
body potential used by us is, however, special and can be evaluated without keeping track 
of triplets. We illustrate below how this may be done for our system p]3| , pT| and derive an 
explicit expression for the energy. 

The three body part of the energy (see Eq. (1) ) is given by, 

^3 = ^1] $3(rij,rjfc,rfci) 

= J2 2(sin2^,cos2^,-4sin^^,cos^e,)^/,fc (12) 

Now define Xij = Xij/vij and yij = Vij/rij so that sin^j = Xikjjij — Xijijik and cos^j = 
XikVij + Xijfjik. Using the above definitions and the quantities. 





~2 ~2 f 


















~2 ~2 ( ~2 
- ^ijVijy^ij 


- 






- yij)fij 



(13) 



We get, ^3 = V^sEi^i with. 



5, = 4[G,(1)F,-4G,(1)2-G,(2)2] - 16 x 

{Gi{3)Fi + 32 Gi{3f + 2 0,(4^ + ^,(1)^ - 

16 G',(3)G,(1) - 4 G,{5)G,{2) + 16 ^,(5)^} (14) 

and Gi{n) = Y,j^i9ij{n) and = Y.j^ifij- The three body forces can be got by taking 
derivatives of which can be cast into a similar form. 
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FIGURES 




FIG. 1. Plot of the second order direct correlation function C^'^\q) vs. wavenumber q for a 

supercooled hard disk liquid in two-dimensions. The lines mark the lengths of the reciprocal lattice 
vectors for the triangular (top) and the square (bottom) lattices scaled so that the smallest RLV 
corresponds to the first peak in C^'^\q). 
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FIG. 2. Definition of angles and distances used in the 3-body potential. 
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FIG. 3. Zero temperature phase diagram in the V3 — p plane. The regions where the square 
and the triangular phases are stable are labelled. 
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(i) / (ii) 

• ; ^ b 

e 

FIG. 4. Two equivalent ways of obtaining a triangular lattice (filled circles) from a square (i) 
and (ii). One can either (i) shear the original lattice by an angle 9 or (ii) rotate the original lattice 
by 45° and then stretch it along one of the axes and compress it along the other so that b/a > 1. 
For a square lattice 9 = 0,b/a = 1 and 9 = 15°, b/a = s/S for the ideal triangular lattice. In terms 
ofthe shear strain £3 the corresponding numbers are and .27812. 
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FIG. 5. Energy difference per particle {AE/N) between the square and rhombic lattices as a 
function of the strain order parameter e (see text). Of the three minima shown, the one at e = 
corresponds to the square phase and the two degenerate minima at e = iteo corresponds to two 
different orientations of the triangular phase. The curves are for V3 = 2.0 (top), 1.5 and 1.0. Note 
the first order transition from the square to the triangular phase as V3 is reduced. 
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FIG. 6. The pressure P = a\ and the "effective" shear stvessdE/de = (73 + Pe as a function 
e for p = 1.1 and V3 = 1.5. 
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FIG. 7. The T"^ order clastic moduli Cn (bulk), C22 and C33 + P (shear) as a function of e. 
Note that for the triangular lattice C22 = (^33 + -P as required by symmetry. The density p = 1.1 
and V-i = 1.5 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

e 



FIG. 8. The per particle Helmholtz free energy difference (AF/N) between the square ((e) = 0) 
and rhombic lattices as a function of (e) for various temperatures T = .1 (top), .5 and .1. Note the 
first order phase transition from square to triangular with increasing temperature. 
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FIG. 9. Phase diagram in the T — V3 plane for p = 1.05 (a) and 1.1 (b). For large V3 the 
square phase is stable while the triangular phase is stable for smaller values of V3. The points 
are results from our molecular dynamics simulations in the NVT ensemble with 2499 particles. 
Starting from an initial ideal square lattice the system was equilibrated for upto 60000 steps and 
the final structure noted (□ for square and A for triangular) for various values of T and V3. The 
solid line is the phase boundary resulting from the cell model approximation (see text) and the 
dashed line is the metastability limit for the square phase from the same theory. 
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FIG. 10. The energy difference AE/N for two values of the anisotropy parameter a = and 
1 for V3 = 1 and p = 1.1. Note that the energy minima for e ^ shifts to lower values of e as a 
increases, (inset) The meanings of the angles 6ij and (p used in the text. While 6ij is the angle of 
the position vector between the molecules i and j measured with respect to the crystal axis {01} 
in the reference square lattice, cf) is the orientation of the basis molecule measured with respect to 
the same axis. 
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FIG. 11. (left)The zero temperature phase diagram in the Vs — a plane for p = 1.1. The 
dashed hne marks the Hmit of metastabihty of the square phase. Note that for V3 = reducing 
a produces a second order transition with a tricritical point at ate = 2.24. The inset shows the 
jump in the order parameter Ae accross the square -rhombus phase boundary as a function of the 
anisotropy a. (right)The phase diagram in the V3 — T plane for a = 1 and p = 1.1. The dashed 
line marks the limit of metastabihty of the square phase. Note that in real systems a depends on 
T so that, in general, any quench traverses a trajectory in the parameter space T — a — V3 and the 
first-order line may end in a non-zero temperature tricritical point. 
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